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Abstract 

Despite the large number of studies available on nicotinic acetylcholine receptors, a complete account of the mechanistic 
aspects of their gating transition In response to llgand binding still remains elusive. As a first step toward dissecting the 
transition mechanism by accelerated sampling techniques, we study the ligand-induced conformational changes of the 
acetylcholine binding protein (AChBP), a widely accepted model for the full receptor extracellular domain. Using unbiased 
Molecular Dynamics (IVID) and Temperature Accelerated Molecular Dynamics (TAMD) simulations we investigate the AChBP 
transition between the apo and the agonist-bound state. In long standard MD simulations, both conformations of the native 
protein are stable, while the agonist-bound structure evolves toward the apo one if the orientation of few key sidechalns In 
the orthosteric cavity is modified. Conversely, TAMD simulations initiated from the native conformations are able to produce 
the spontaneous transition. With respect to the modified conformations, TAMD accelerates the transition by at least a factor 
10. The analysis of some specific residue-residue interactions points out that the transition mechanism Is based on the 
disruption/formation of few key hydrogen bonds. Finally, while early events of ligand dissociation are observed already In 
standard MD, TAMD accelerates the llgand detachment and, at the highest TAMD effective temperature, It is able to 
produce a complete dissociation path in one AChBP subunlt. 
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introduction 

The transmembrane nicotinic acetylcholine receptors 
(nACliRs), belonging to the so-called 'Cys-loop' super-family of 
ligand-gated ion channels (LGICs) [1,2], are involved in a variety 
of biological processes [3-5] and have been implicated in the onset 
of Alzheimer's disease [4,6] and nicotine addiction [7]. The 
nAChR channel pore, located in the transmembrane domain 
(TMD), opens following the binding of agonist ligands in the 
orthosteric site located in the extracellular domain (ECD) of the 
protein. Conversely, the channel is mainly closed in the resting 
state, when either no ligand is present or an antagonist is bound, 
and in the desensitized, agonist-bound state. The main endoge- 
nous agonist for iiAChRs is acetylcholine (ACh). Nicotine also 
binds nAChR as an agonist while lobeline (Fig. SI) has a full or a 
partial agonist efTect [8,9]. 

Currently, only limited information is available on the atomic 
structure of nAChRs. Low resolution structures of Torpedo 
acetylcholine receptor with closed [10] and open [11] pore were 
obtained from electron microscopy data. Structures of distant 
prokaryotic homologues of nAChRs, present in Gheobacter violaceus 
and in Erwinia chrysamthemi, GLIC and ELIC, were solved during 



the last years by high-resolution crystallography [12-17,17-19]. 
Interestingly, the GLIC structures put in evidence a binding site 
for anesthetics in the TMD [20], ketamine binding to the 
orthosteric site [21] as well as a new semi-closed state of the 
channel [22]. The ELIC structures showed the binding of 
anesthetics both in the TMD and in the ECD [23] and the 
binding of acetylcholine to the ECD [2 1] . However, despite these 
new pieces of information, the nAChR gating conformational 
transition is not yet fuUy understood. Valuable insight on the 
ligand binding mechanism came from studies of a water-soluble 
homologue of the nAChRs ECD, the pentameric acetylcholine 
binding protein (AChBP). Indeed, AChBPs have been crystallized 
bound to different ligands displaying agonist or antagonist effects 
on iiAChRs [24,25]. The structures of unliganded (apo) [26,27] 
and hganded (holo) state of AChBP [27-40] revealed how the 
ligands bind to the orthosteric pocket. This pocket, present in each 
of the subunits, is lined by the so-called loop C and extends at the 
interface between subunits (Fig. 1). The AChBP structures 
revealed the influence of the ligand type on the degree of C- 
loop closure against the protein core [32], as the loop arrange- 
ments cluster into three groups: i) the antagonist-bound "open" 
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Principal subunit Complementary subunit 



Figure 1. The structure of AChBP. The five homomeric subunits in AChBP are labeled SI, S2, S3, S4, S5, in a clockwise direction as viewed from 
the apical side and colored black, red, green, blue, yellow respectively. A) View from the apical side of the AChBP protein, shown as a cartoon model; 
B) Side view; C) The AChBP subunit interface bound to the ligand lobeline, represented in sticks. The relevant loops on both the principal and 
complementary subunit are indicated. 
doi:1 0.1 371 /journal.pone.0088555.g001 



conformation; ii) the apo "intermediate" conformations; iii) the 
agonist-bound "closed" conformations. It is worth noting however 
that a recent AChBP structure [41] questions this model, in the 
case of the small antagonist ligand dihydro-beta-erythroidine 
(DHjSE). 

Although the ECD/TMD interface region [31], involved in 
discriminating between agonist versus antagonist binding and in 
the long-range communication between the extracellular binding 
site and the transmembrane pore, is not present in AChBP, similar 
structural patterns and pharmacological responses have made 
AChBP a valuable proxy for investigating ligand recognition by 
nAChRs. 

Due to the lack of high-resolution structural information on the 
nAChRs, computational methods have been extensively used to 
probe the conformational transition. Nanosecond time scale 
Molecular Dynamics (MD) simulations have been performed on 
homology models of nAChRs [42-44] . Along with normal mode 
analysis, performed by using an elastic-network representation for 
the protein [45,46] or on a fuU atomistic protein model [44,47,48], 
they provided first useful insights into the gating mechanism. 
However, because the binding-to-gating process in nAChRs takes 
place in the millisecond timescale in physiological conditions [49], 
standard MD approaches are not useful. Targeted MD has been 
used, in which the C-loops have been forced to move from the 
outward to the inward, agonist bound, position [50,51]. 10 ns MD 
simulations have been performed on AChBP in complex with 
acetylcholine [52], providing a picture on the ligand interaction 
with binding site residues; with nicotine and carbamylcholine [53], 
focusing mostly on the dynamics of residues and water in the 
binding pocket rather than on large global motions, and with 
partial agonists [54]. Spontaneous, although preliminary, confor- 
mational changes have been observed in unbiased MD simulation 
of a homology model of the (3(7 nACliR ECD both in the apo state 
and bound to an antagonist toxin [55]. Steered MD was used to 
study the unbinding of nicotine from the AChBP protein [56]; the 



unbinding of acetylcholine from a homology model of the ECD of 
the nicotinic receptor has been investigated, by using as reaction 
coordinate the distance between the ligand and the binding site 
[57]. 

In the present work we use full-atom standard MD and tiie 
enhanced sampling method temperature-accelerated MD 
(TAMD) [58] to study the AChBP holo-to-apo transition. In 
particular, we investigate the conformational changes involved in 
the opening and closing of the C-loops. In TAMD, a set of extra 
variables are introduced, coupled to the original system via 
collective variables (CVs). The original and new variables are then 
evolved concurrentiy in condition of effective adiabatic separation 
and at two different temperatures, the physiological temperature 
for the physical system and a higher artificial temperature for the 
new variables. In this way, the system navigates the free energy 
landscape associated to the extra variables overcoming barriers 
that are even much higher than the energy at the physiological 
temperature. TAMD can be used to reconstruct the free energy 
landscape from direct sampling via reweighting [58,59], or 
combined with a mean-force interpolation method [60]. TAMD 
has already been applied to a variety of rare events sampling 
studies [61-67] and has proved to be particularly useful in protein 
conformational searches [68-75]. Here, we use it to search for 
AChBP conformations by exploring the space of a set of suitable 
variables. 

The conformational transition between the lobeline-bound 
(holo) and the apo states of Aplysia califomica AChBP is investigated. 
First, unbiased MD simulations on the hundred of nanoseconds 
time scale allowed to characterize the relative mobilities of the 
different AChBP regions and to give an insight on which CVs are 
suitable to describe the transition between the holo and apo 
conformations. MD simulations starting from ad hoc modified holo 
conformations put in evidence the influence of few key residues on 
the nature of the metastabUity, or in other words, let us to locate 
the transition bottlenecks at molecular level. Then, TAMD 
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simulations were performed using as CVs the centers of mass of 
the C-loop and the cys-loop, to accelerate sensibly the holo-to-apo 
transition of AChBP. 

Results 

MD Simulations 

MD simulations of AChBP were carried out for several systems 
in explicit water, each 150 ns in length: the apo pentamer (P), the 
lobeline-bound pentamer in the presence (Pl+L) or in the absence 
(PI) of lobeline. Details on the starting structures and simulation 
methods are given in Section "Materials and Methods" and in 
Text SI. 

As first, an analysis of the protein secondary structure has been 
performed along the trajectories. The content in a. helices and /? 
strands, calculated using STRIDE [76], mostly superimpose with 
the corresponding contents in the X-ray crystaUographic structures 
(data not shown) and is quite similar in all trajectories and subunits 
except for the two small 3-10 heUx (residues 61-63) and o(-helix 
(residues 67-70). While the shortest disappears along the 
simulations, the residue stretch at 67-70 still shows a significant 
percentage of 3-10 helix content along the whole trajectory in 
each simulation done. We attribute these secondary structure 
changes to the absence of crystal packing. 

To further assess the stability of the model systems along the 
trajectories, the Root Mean Square Deviation (RMSD) of the 
individual subunits with respect to the starting conformations after 
equilibration was calculated after removing the roto-translational 
motions of each individual subunit (Fig. S2). All profiles show a 
plateau after 20 ns, at about 1.2 A for Pl+L and 1.7 A for P and 
PI. In Pl+L, the RMSD curves of all subunits are superimposed 
within 0.5 A, whereas they are superimposed within 1.0 A and 1.5 
A in P 1 and P, respectively: this means that the absence of ligand 
cdlows larger conformational changes and more variability among 
subunits. This is in agreement with results from previous 
simulations of AChBP [53]. 

The atomic Root Mean Square Fluctuations (RMSFs) profiles 
(Fig. 2) are similar for the three trajectories Pl+L, PI and P, with 
the exception of few regions: the C-loop, the F-loop, the al-/?l 
loop (residues 9-27) connecting the helix al and the strand /?1, 
and the ^1-^1 loop. With respect to P, in Pl+L die RMSFs 
averaged among subunits are smaller in these regions. In 
particular, the loop C displays the largest decrease, followed by 
the loop F. Indeed, because of the interaction with lobeline the 
loop C is locked in a specific conformation with limited 
fluctuations [27]. On the other hand, the fluctuations in cys-loop 
region displays quite similar values in all simulations. Despite the 
overall similar RMSFs profiles in trajectories Pl+L, PI and P, the 
differences between P 1 and P show that removing ligand does not 
produce the same RMSFs profile obserx^ed for P in the time 
window considered. The large difference in the C-loop fluctuations 
is due to the opening of C-loop in P (the distance defined in 

Section "Trajector)- Analysis", is 15.1 A, a\Tragcd o\'er the 
trajectory and over the five subunits) to a larger exti;nt than in P1+ 
L (averaged </g'™=13.5 A) and PI (averaged rfg'™=13.8 A). 
Differences in opening correspond to those observed in the starting 
X-ray crystaUographic structures: 2BYS (rf""'" = 12.7 A, averaged 
over the five subunits) and 2W8E (averaged i^g'™ = 14.7 A). 

The cys-loop does not display significant displacement with 
respect to the remaining part of the pentamer, as the distances 
and d°f^' defined in Materials and Methods (Section 
"Trajectory analysis") display variations smaller than their 
corresponding standard deviations (data not shown). This is most 



likely related to the structure of AChBP, which encompasses only 
the ECD of nicotinic receptors. 

Distributions of were calculated along the trajectories, by 
joining together data from all five subunits. Results are shown in 
Fig. 3. The Pl+L distribution is unimodal, centered on the value 
13.4 A, whereas the P distribution is bimodal with two peaks 
located at 13.2 and 15.6 A (Fig. 3, first row). The PI distribution is 
also uni-modal, although broader than Pl+L and displaying a 
small shoulder at 12.4 A. The distributions calculated along the 
three time intervals: 0-50, 50-100 and 100-150 ns (Fig. 3, second 
row) are almost superimposed for Pl+L and P whereas they 
display more heterogeneity for PI. This is the sign of a reached 
stable condition of the Pl+L and P trajectories. PI is not yet in 
stable conditions and this is probably due to the destabilization 
produced by removing lobeline. This destabihzation is not 
overcome in 150 ns, as is also visible in RMSD and RMSFs plots 
(see Fig. S2 and Fig. 2). 

Starting from biased conformations of AChBP (confl and conf2, 
described in Text SI, section "New AChBP models") removes this 
destal)ilizati<)ii and makes the c?™"'" distribution evolve along 
100 ns toward the bimodal distribution observed in P (Fig. 3, third 
and fourth rows). The drift is even more pronounced in Plconf2, 
where the LYS143 sidechain was modified. At variance, in the 
presence of lobeline (Plconfl+L, Plconf2+L, see Table 1), the 
fjintra (Jijtj-ijjufion Stays unimodal as in th(- simulation Pl+L. The 
presence of lobeline is thus able to counterbalance the introduced 
perturbation, even if the orientation of LYS143 is changed. 

To summarize, no major change in the average pentamer 
architecture is observed along the MD trajectories of Pl+L, PI 
and P. Larger conformational drifts and heterogeneity in atomic 
fluctuations appear in the absence of lobeline and in the apo 
system P. The distribution of C-loop openings in Pl+L/Pl on one 
side, and in P on the other side, are respectively unimodal and 
bimodal, revealing the existence of two metastable states which 
cannot easily inter-convert in the MD timescale here considered. 
Biasing the initial sidechain conformations of PI induces the C- 
loop opening, in particular if LYS143 is perturbed. 

TAMD Simulations 

The data discussed so far confirm that the closed to open 
transition of the C-loop occurs on very long time scales. To 
accelerate this transition, three sets of TAMD simulations at four 
different effective temperatures were carried out using systems all 
with closed C-loop conformations (see Table 1). The first two sets 
are for unliganded AChBP and were started from the conforma- 
tions corresponding to the 40 ns snapshot of the PI and Pl+L free 
MD simulations, removing the ligand in the second case. We 
indicate them as (PI), and (PI -40 ns);, where / indicates the 
effective energies equal to 3, 5, 7 and 10 kcal/mol. The third 
set is for liganded AChBP, starting from the 40 ns snapshot of the 
Pl+L free MD and is indicated as (Pl+L),. 

As a first check on the TAMD simulations we calculated the 
protein secondary structure contc-nt along the trajc-ctorics. 
Similarly to what was observed along the standard MD, the a 
helix and J? strand content along all TAMD trajectories mostiy 
superimpose with the corresponding contents in X-ray crystaUo- 
graphic structures (data not shown), with the exception of the two 
small 3-10 helix (residues 61-63) and a-helix (residues 67-70), as 
already observed along the MD trajectories. This points out that 
during TAMD the protein secondary structure is not altered, with 
respect to the MD. 

In order to check if TAMD has been successful in driving the C- 
loop opening, the t/"'™ distributions along the TAMD trajectories 
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Figure 2. Root Mean Square Fluctuations of Ca atoms. Root Mean Square Fluctuations (in A) of Ca carbon atoms in the individual subunits 

/V"' (rc,(l;)-<rr,»- 

with respect to the average structure calculated as RMSF(Ca)= \J ' — , where tcStj) ^^i^ position of the Q atom at yth time step, N, 

is the total number of steps and (r^) = jf-^Yl'jLx fc.C'j). along the the 50-150 ns time interval (i.e. where the RMSDs are stabilized, see Fig. S2). Upper 
panel: P; central panel: PI; lower panel: P1+L The curves are colored according to the scheme in Fig. 1. The line in magenta is the average over the 
five subunits. 

doi:10.1371/journal.pone.0088555.g002 



were calculated and compared with the bimodal distribution 
obtained in the standard P simulation (see Fig. 3). Results are 
shown in Fig. 4 for (Pl)„ (Pl-40ns),- and (Pl+L),- from top to 
bottom and ;6^' equal to 3, 5 and 7 kcal/mol. Results for 
equal 10 kcal/mol are shown in Fig. S3. 

It can be observed that in the isolated AChBP system at 
equal to 3 kcal/mol, (Pl)3, the profile of the rfg'™ distribution 
approaches the one for P already in a time scale of few 
nanoseconds (Fig. 4, left upper panel, red bold line). At later 
times, the distribution becomes monomodal, but peaked at larger 
values of d'^^'"' (green and cyan bold lines). Results from TAMD 
started from the PI -40ns structure also show a bimodal 
distribution in the interval 0-4 ns (Fig. 4, left middle panel, black 
bold line); a bimodal distribution even more similar to the one in P 
is achieved in the interval 12-16 ns (magenta bold line). At later 
times, the distribution shifts to higher values of d'^"" (green bold 
line). Thus, TAMD accelerates the opening of the C-loops in the 
unliganded AChBP by at least a factor of 10 with respect to the 
standard MD simulations, already at 3 kcal/mol. 

Conversely, as shown in Fig. 4, left lower panel, the presence of 
the lobeline prevents this acceleration at the smallest value of ^ ' , 



suggesting that a higher effective energy is needed to open the C- 
loops in the presence of the ligands; indeed, the d'^"" distribution 
in the presence of lobeline starts to have a bimodal shape (although 
stiU with higher percentage of closed C-loops) at 5 kcal/mol, on 
the 0-4 ns time scale (Fig. 4, central lower panel, black bold line). 
A secondary small peak at d'^i*'" about 1 7 A appears later, in the 
time interval 4—8 ns (red bold line). Increasing the effective energy 
up to 7 kcal/mol produced an unimodal "open" distribution 
already before 4 ns (Fig. 4, right upper and middle panels); at 
variance, in the presence of lobeline, a bimodal distribution is 
achieved over 0-4 ns (Fig. 4, right lower panel, black curve), with 
roughly the same percentage of closed and open monomers. Then 
the distribution switches to unimodal shape over 4—8 ns, with a 
peak corresponding to an intermediate conformation. A bimodal 
distribution is observed again, over the last 16-20 ns of TAMD 
simulation, this time corresponding to mostly open subunits. 

An almost wide unimodal distribution is observed over 0-2 ns at 
10 kcal/mol in the unliganded case (see Fig. S3, upper panel); in 
the presence of lobeline, bimodal distribution is achieved over the 
first 2 ns (Fig. S3, lower panel). Note, however, that, already at 7 
kcal/mol, the rf^'"'" distributions quickly evolve toward a profile 



PLOS ONE I www.plosone.org 



4 



February 2014 | Volume 9 | Issue 2 | e88555 



Simulation of AChBP 



,1 1— -TIT 

Piconfl 




10-100ns 







P1conf2 



10-100ns 




lO -T 




PIconC+L 



10-100ns 




10 12 14 16 18 20 10 12 14 16 18 20 10 12 14 16 18 20 10 12 14 16 18 20 



5"- 



-I 1 1 1 r 

P1conf1 



10-50ns 
50-100ns 




P 1 1 1 1 r 

P1oonf2 



10-50ns 
SO-IOOns 



1 1 r- 

Piconfl *L 



50- 100ns 



T 1 1 1 1 1: 

PIconC+L 




10-50ns 
SO- 100ns : 



10 12 14 16 18 20 "^10 12 14 16 18 20 "^10 12 14 16 18 20 "^10 12 14 16 18 20 

dC^lA) dC'""(A') AC"'"(k) dC"""(A) 



_1 1 


1 1 1. 

P 




0-15GrB 






1 1 1 



10 12 14 16 18 



r T 1 


1 — 1 1 
p 


r 




0-50ns 






50-100ns 






100-1 sons 






1 " 


1 



Its T? 



"14 W 
dC'"'"(A") 



T8 20 




_i 1 


1 1 

PI 


l_ 




^— 0-150n3 




/ 


V 




: 1 1 II 



20 10 
in 



T? T6 — W 

dC"**(A) 



12 



14 16 18 20 



PI 



0-50ns 

50-100ns 

tOO-ISOns 




1? T6 18 20 

dC'"""(X) 



Figure 3. Distribution of rf™'"' values. First and third rows: distribution of rf""" values (A) calculated over all the MD trajectories; second and 
fourth rows: distribution of J""" values (A) calculated in different time intervals (see legends), d'l^"" values are collected from the five subunits. 
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much more difTerent than the one observed for P, with d'^"" values 
up to 36 A. The presence of lobeline is no more an obstacle to the 
C-loop opening, but TAMD seems to evolve the system to states 
distinct from the natural apo state P. 

Hydrogen Bonds between AChBP Residues 

Let us discuss now the role of specific residue-residue 
interactions in the stabilization of the different C-loop conforma- 
tions. The comparison of the X-ray crystaUographic structures 
2W8E (apo AChBP) and 2BYS (holo AChBP) points out that the 
C-loop closed conformation in the presence of lobeline is stabilized 
by a network of hydrogen bonds, among which is the one between 
LYS143 in /? strand 7 and TYR188 in the C-loop (see Fig. S4F). In 
Fig. 5 we report the persistence of this hydrogen bond along our 
standard MD simulations, defined as the percentage of time along 



a fuU trajectory that the bond is formed. For each simulation, data 
are plotted for the five different protein subunits. As shown in 
Fig. 5A, the LYS143/TYR188 bond is very rarely formed in P, 
while it is more formed in most of the subunits in the presence of 
lobeline. Modifying the initial PI protein structure makes the 
hydrogen bond percentage evolve toward the one in the native apo 
P state, particularly for PlconK. These structures relaxes plausibly 
toward the target state, i.e. the unbound one, on the time scale of 
the unbiased MD simulation, pointing out where the transition 
bottlenecks could be located, at molecular level. Along the TAMD 
simulations started from the PI and PI -40 ns initial conditions 
(Fig. 5B) the behavior is much more similar to the one observed in 
the target P protein with respect to the one along the standard MD 
PI. TAMD produces the evolution toward P spontaneously 
already at equal to 3 kcal/mol, and even more at equal 
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Table 1. Simulations of the AChBP pentamer. 



Name 


Starting point 


Components 


Type 


^ ' (Itcal/mol) 


Duration (ns) 


Pl+L 


PDB 2BYS 


AChBP, LOB 


MD 


- 


150 


PI 


PDB 2BYS 


AChBP 


MD 


_ 


150 


P 


PDB 2W8E 


AChBP 


MD 


- 


150 


Plconfl+L 


PDB 2BYS" 


AChBP, LOB 


MD 


- 


100 


Plconf2+L 


PDB 2BYS'' 


AChBP, LOB 


MD 


_ 


100 


Plconfl 


PDB 2BYS" 


AChBP 


MD 


_ 


100 


Plconf2 


PDB 2BYS'' 


AChBP 


MD 


_ 


100 


(P1)3 


PI at 40 ns 


AChBP 


TAMD 


3 


30 


(Pl)5 


PI at 40 ns 


AChBP 


TAMD 


5 


30 


(Pl)7 


PI at 40 ns 


AChBP 


TAMD 


7 


30 


(Pl)io 


PI at 40 ns 


AChBP 


TAMD 


10 


10 


(Pl+U, 


Pl+L at 40 ns 


AChBP.LOB 


TAMD 


3 


20 


(Pl+L)5 


Pl+L at 40 ns 


AChBP, LOB 


TAMD 


5 


30 


(Pl+L), 


Pl+L at 40 ns 


AChBP, LOB 


TAMD 


7 


20 


(Pl+L)io 


Pl+L at 40 ns 


AChBP, LOB 


TAMD 


10 


10 


(Pl-40ns)3 


Pl+L at 40 ns'' 


AChBP 


TAMD 


3 


30 


(Pl-40ns)5 


Pl+L at 40 ns'' 


AChBP 


TAMD 


5 


20 


(Pl^Ons)? 


Pl+L at 40 ns' 


AChBP 


TAMD 


7 


20 



"New models (see Section: "New AChBP models" in SM); 
''as before, plus disrupting the LYS143/TYR188 bond. 

'The ligands were removed after 40 ns of the Pl+L unbiased simulation. LOB: abbreviation for Lobeline. 
doi:l 0.1 371 /journal.pone.0088555.t001 



-kgT where T is the artificial temperature. 



to 5 kcal/ mol. We would like to stress that TAMD is not used here 
to validate the method against the simulation results on the ad hoc 
structures, as it could appear, but to explore the free energy 
underlying the C-loop transition, with the aim to overcome the 
known bottlenecks. A posteriori, the results obtained remark the 
usefulness and powerfulness of the TAMD method per se. 

The relation between the LYS143-TYR188 hydrogen bond and 
the degree of C-loop closure is better evidenced by plotting the two 
dimensional histograms of d'^"" and the distance between LYS 1 43 
and TYR188 donor/acceptor atoms. In Fig. 6 we show such maps 
calculated over aU unbiased MD trajectories (panels A-C). When 
the bond is formed or partially formed, as in PI and Pl+L (Fig. 6, 
panels B and C), the C-loop is in the closed state; in 
correspondence to larger LYS143-TYR188 distances as in P 
(panel A), the C-loop is mainly in the open form. 

The strength of the LYS143-TYR188 bond and the possibility 
that its formation may contribute to the C-loop closure was 
already investigated for the apo AChBP in [50], where the 
potential of mean force for the interaction of the two above 
mentioned residues was calculated using umbrella sampling on the 
LYS143NZ-TYR1880H distance. In that work, the bound-like 
conformation was stabilized by about 2.0 kcal/mol at a distance of 
3A. The strength of the hydrogen bond was found somewhat 
weak; according to the authors, this could arise from the flexibility 
of the C-loop and/or concurrent interactions of both residues with 
water molecules. 

At variance with monodimensional PMF calculation, we study 
the transition in the fuU protein, i.e. we accelerate all C-loops 
together. According to our results, the formation and persistence 
in time of the bond between LYS 143 and TYR 1 88 is strictly anti- 
correlated to the presence of the lobeline. The hydrogen bond is 
weak and transiently formed in Pl+L, as distance values larger 



than 3A are observed (see Fig. 6, panel C), while the distribution in 
PI show a larger population in the range 2-4 A (Fig. 6, panel B). 
The weakness in Pl+L is brought about by the transient hydrogen 
bond that the lobeline could estabhsh with TYR 188 (see below). 
Therefore, in AChBP bound to lobeline, the hydrogen bond 
between LYS 143 and TYR 188 does not have a dominant 
structural role in maintaining the C-loop closed. 

The relation between the LYS143/TYR188 donor/acceptor 
distance and J™"''' was followed also along the TAMD trajectories. 
Results are shown in Fig. 6, for the TAMD starting from PI and 
PI -40ns at equal to 3 and 5 kcal/mol (panels D to G), 

confirming i) how TAMD was able to efiiciendy push the PI 
protein structure out of the basin in which it was trapped along the 
150ns MD trajectory (see Fig. 6, panel B), already at equal to 
3 kcal/mol (see Fig. 6, panel D); ii) the higher energy cost for 
opening the C-loops and disrupting the LYS143/TYR188 bond 
when the lobeline is present in the binding pocket, as pointed out 
by the comparison of the maps in panels H, I. 

Another hydrogen bond that has been long discussed for 
AChBP is the one between SER146 and TYR93 (see Fig. S4A). It 
has been suggested that TYR93 has a role as "gatekeeper" in 
AChBP [77]: depending on its position, it makes the lobeline 
pocket inaccessible (as in 2W8E) or open (as in 2BYS). 
Furthermore, experimental results on fuU nicotinic receptor [78] 
suggest that mutating TYR93 (in particular aY93W) causes a 
change in the structure of the binding pocket that compromises 
both the mobility of the hgand into and out of its docking site, and 
the speed of channel opening. 

In Fig. 7 we report the persistence in time, along our 
simulations, of the SER146/TYR93 bond: in P the bond is 
always present (Fig. 7A), while in Pl+L it is never formed; this is 
due to the presence of the lobeline in the binding pocket. As it wiU 
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be shown below, the carbonyl group of SER146 interacts with tlie 
lobeline hydroxyl group along the entire Pl+L trajectory, in four 
out of five subunits. In the unbiased PI trajectory, an intermediate 
behavior between P and Pl+L is observed. Once again, biasing the 
initial PI protein structure makes the hydrogen bond percentage 



to evolve toward the one in the natural P conformation, 
particularly for PlconC 

Along the TAMD simulations started from the PI and PI -40ns 
initial conditions (Fig. 7B) the persistence in time becomes much 
more similar to the one observed in the target P protein with 
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Figure 5. Hydrogen bond analysis. Percentages of hydrogen bond formation between LYS143 and TYR188 in each subunit along MD (A) and 
TAMD (B) trajectories. The hydrogen bond is considered to be formed when the minimum distance among all possible donor/acceptor distances, 
between the two residues, is less than 2.5 A. 
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Figure 6. Two dimensional Kiistograms. Two dimensional histograms n{d'^'"', LYS143/TYR188 closest donor/acceptor distance). A: simulation P; 
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3 kcal/mol and 5 kcal/mol; H and I: TAMD from P1+L initial condition at 3 kcal/mol and 5 kcal/mol. 
doi:10.1371/journal.pone.0088555.g006 



respect to the one along the standard MD PI. TAMD produces 
the evolution toward P spontaneously already at ;8^' equal to 3 
kcal/mol. In analogy with the LYS143/TYR188 bond, the two 
dimensional histograms shown in Fig. 8 point out a relation 
between close/ open C-loops and the formation/ disruption of the 
SER146/TYR93 bond. It is clearly shown how TAMD was 
efficient in driving the system toward the natural P state already at 
j8 ' 3 kcal/mol; this was better achieved by starting from the Pl- 
40ns initial condition. 

The features of the SER146/TYR93 bond in the TAMD 
exploration are further discussed in Text SI (see section " T/^MD 
exploration and the role of SER146/TYR93", Fig. S8, Fig. S9 
and Fig. SIO). 

AChBP-lobeline Interaction 

We now turn to discuss the lobeline /AChBP interactions and 
conformations along our MD and TAMD simulations of the 



complex systems. At first, we monitored [27]: (i) four hydrogen 
bonds connecting lobeline carbonyl and the hydrogen atom Hfil 
of TRP147, lobeline amide and carbonyl of TRP147, lobeline 
hydroxyl and carbonyl of SER146, and lobeline hydroxyl and 
carbonyl of TRPI47; (ii) three van der Waals interactions: between 
lobeline piperidine ring and TRP147 indole, between lobeline 
methyl and TYR188 aromatic ring, and between lobeline methyl 
and TYR195 aromatic ring. The hydrogen bonds are monitored 
through the distances between hydrogen donor and acceptor 
atoms. The van der Waals interactions are monitored through the 
distances between the centers of mass of the two atom groups. 
Among the hydrogen bonds observed in the crystaUographic 
structure 2BYS [27], only one, between the lobeline hydroxyl 
group (LOB-OH) and the SER146 carbonyl (SER146-CO) is 
formed during the MD simulations (Table SI). Between the rings of 
TYR188 and TYR195, the lobeline methyl (LOB-CH3) seems to 
prefer to interact with TYR195 ring, as the distances LOB-CH3/ 
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TYR195-ring are smaller than the ones observed for LOB-CH3/ 
TYR188-ring (Table SI). This marked difTerence from the X-ray 
crystallographic structure can be explained by the fact that only the 
global electronic envelop of the ligand was visible in the Xray 
crystallogTaphic structure 2BYS, thus hampering a precise defini- 
tion of the ligand conformation. Among the protein-lobeline 
distances monitored on the trajectory Pl+L, in the subunit P3, 



the distance between the lobehne hydroxyl (LOB-OH) and the 
carbonyl groups of SER146 and TRP147, displays a strong 
transition around 64 ns (Fig. 9A, green and red curve, respectively). 
This transition is accompanied by the formation of a hydrogen bond 
between the lobeline amide (LOB-NH) and the carbonyl group of 
TRP147 (Fig. 9A, black curve). An examination of frames extracted 
in the 65-100 ns range of Pl+L reveals that one of the aromatic 
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Figure 8. Two dimensional histograms. Two dimensional histograms n(d'^"", SER146/TYR93 closest donor/acceptor distance). A: MD simulation 
Pl+L; B: PI; C: P; D: TAMD simulation (Pl+L)^; E: (PI),; F: (P1-40ns)3. 
doi:1 0.1 371/journal.pone.0088555.g008 
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rings of lobeline changed its rotameric state around 64 ns, and the 
LOB-OH tlien pointed toward the complementary subunit, as 
represented in Fig. 9D, establishing transient hydrogen bonds with 
the hydroxyls of the TYR55 (Fig. 9A, purple curve), TYR93 
(Fig. 9B, black curve) and TYR188 (Fig. 9B, green curve). This 
movement of lobeline is an isolated event, which might describe the 
initial steps of lobeline dissociation. 

The ligand lobeline has been shown to have a partial agonist 
effect [79], on a7 [8] and a4/?2 [9] nAChRs. Although it is not 
clear what partial agonism means at the molecular level [80] , the 
observation of the lobeline transition would support a model 
where a partial agonist is an agonist ligand that interacts weakly 
with the receptor and is thus more prone to dissociation. On the 
other hand, it has been observed that other partial agonists co- 
crystallized with AChBP cause an incomplete closure of loop C 
against the binding site [32]; at variance, in 2BYS the C-loop is 
capped around the lobeline, similar to what observed with fuU 
agonists. This has been also observed in the recent Xray 
crystallographic structure of AChBP from Capitelh teleta, bound 



to lobeline [81]. Our results (see the c/™'™ histogram in Pl+L in 
Fig. 3) show that the lobeline-induced tight closing is kept during 
the simulations. The lobeline OH transition would reconcile these 
two points of view as the lobeline weakened interaction would arise 
from isolated dissociation events rather than from different, 
induced, conformations of C-loop. 

In the biased MD simulations Plconfl+L and Plconf2+L, the 
lobeline transition was not observed. This could be explained by 
considering that the starting points were produced by modifying, 
among others, the orientation of GLN186, which strongly 
influences the position of TYR188 (data not shown). 

The lobeline transition is however observed in TAMD 
trajectories (Table S2), and significantly accelerated in at least 
one subunit (see Fig. 10 for a representative subunit in each 
TAMD trajectory). In particular, it is observed progressively in 
more than one subunit as the effective temperature increases. At 
equal to 3 kcal/mol it occurs in one subunit only, namely P3 
at about 15 ns. At 5 kcal/mol it occurs in three out of five 
subunits: PI first at about 15 ns, then P4 at about 22.5 ns and 
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Figure 9. The ligand transition. A and B: residues-lobeline donor/acceptor distances (see legends); C: lobeline and protein residues configuration 

before the OH transition; D: after the OH transition. 
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finally in P5 about 27 ns. At equal to 7 kcal/mol, it occurs in 
four out of five subunits: P4 first at 2.5 ns, then P3 at about 5 ns, 
followed by P5 at about 10 ns and finally in PI at about 18 ns. At 
the highest 10 kcal/mol, the event is observed in the five 

subunits at the following time frames along the trajectory: P4 at 
about 0.6 ns, P3 at about 3 ns, P2 at about 3.3 ns, P5 at about 
6.4 ns, and PI at about 6.6 ns. Remarkably, we observe that the 
same chain of events characterizes the lobeline transition in MD 
and TAMD, i.e. we observed the same key residues bonds 
rupture/formation events (see Figs. 9 and 10). 

In the TAMD at ;6^' equal to 10 kcal/mol, the dissociation of 
the lobeline from the binding pocket is complete after few 
nanoseconds (Fig. 11). The ligand exits from the binding site and 
migrates downward to the solvent; its motion is assisted by a 
corresponding downward motion of the C-loop. A similar path 
was obtained for the unbinding of acetylcholine from the ECD of 
a7 nicotinic receptor using Steered MD in [57]. In that work, the 



free energy difference between the bound and unbound states was 
estimated to be about 10-15 kcal/mol, consistent with the effective 
energy adopted by us. At variance with the paths proposed by 
Zhang et al. [57] however, the dissociation path observed in the 
present work was not obtained by forced pulling but occurred 
spontaneously in the TAMD trajectory. A downward motion like 
the one we observed would be not supported by the analysis of 
several crystal structures of AChBP both in the holo and the apo 
state, in which it is shown that the loop C assumes different 
intermediate positions mainly in the direction outwards perpen- 
dicular to the main axis of the protein. We would like to stress that 
the ligand dissociation trajectory observed here is a possible one as 
it is provided by TAMD. Indeed, TAMD trajectories are likely to 
pass close the most probable reaction path [69] although there is 
no explicit requirement for this. In particular, the first step of 
lobeline dissociation is additionally supported as the most plausible 
part of the dissociation because the same event is spontaneously 
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Figure 10. Protein-lobeiine donor/acceptor distances. Protein-lobeline donor/acceptor distances involved in the ligand transition in some 
representative subunits along the TAMD trajectories. First row from top:(P1+L)3 trajectory; second row: (Pl+Us; third row: (Pl+Uy; fourth row: (P1 + 
L)io. Left column: LOB-NH/W147-CO (black), LOB-OH/SI 46-CO (green), LOB-OH/W147-CO (red) and LOB-OH/Y55-OH (purple) hydrogen bond 
distances. Right column: LOB-OH/Y93-OH (black) and LOB-OH/Y188-OH (green) hydrogen bond distances. 
doi:1 0.1 371/journal.pone.0088555.g01 0 
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Figure 11. Snapshots along the TAMD trajectory. Snapshots tal<en at 1 ns, 3 ns, 4 ns, 5 ns, 6 ns and 6.3 ns along the (P1+L)io trajectory for S3- 
S2 interface. S3 is in green; S2 is in red, cartoon representation; lobeline in VDW representation. Lobeline is fully in the solvent at 6.3 ns. 
doi:10.1371/journal.pone.0088555.g011 



observed along the MD simulation of the native lobeline-bound 
complex. 

Discussion 

In the present work, the transition of the Apljda califomica 
AChBP from the hole to apo state was explored using MD and 
TAMD simulations. The transition has been qualitatively 
described using: (i) opening distances of C-loop, (ii) individual 
hydrogen bonds. The TAMD simulations have been shown to 
dramatically accelerate the transition observed in some MD 
simulations only after some ad hoc modifications of the structure. 
Indeed, to get the holo/apo transition of the loop C in unbiased 
MD simulations it was needed to modify the native structure of the 
apo protein (PI), by disrupting the pattern of some key hydrogen 
bonds. Possible first steps of spontaneous lobeline dissociation by 
TAMD have been obseived. At the highest effective temperature, 
a path of complete ligand dissociation has been observed in one 
subunit. 

The analysis of MD simulations has shown that the AChBP 
states bound to lobeline and apo are stable within an interval of 
150 ns. This agrees with the timescale of channel opening in 
nAChRs, which is on the order of 1 ms [49]. AChBP bound to 
lobeline shows along the MD trajectories limited oscillations 
around a well defined mean structure. On the contrary, the apo 
AChBP obtained from 2W8E has a much less defined conforma- 
tional state due to the large C-loop internal mobility, although 
nAChRs in the corresponding apo state display a well-defined 
physiological state with a closed channel. 

The initial change of orientation of few residues (ARG57, 
ARG95, GLN103, GLN184, LYS143) is sufficient to reduce, at 
least in part, the long metastabUity of AChBP in standard MD 
trajectories. These residues have been selected as they display a 
variability of orientations in the AChBP crystal 2BYS. This 
variability in the crystal is probably the sign of an equilibrium of 
AChBP in solution among several conformations. It is interesting 
to see that pushing the conformations of these mobile residues 
toward the ones observed in the absence of lobeline (PDB entry: 



2W8E), induce a whole transition of AChBP toward the apo state, 
as seen by monitoring the C-loop behavior. Thus, the relevance of 
these residues modifications is supported by the plausible effects 
the modifications have on the ACHBP. One should notice that, 
among the modified residues, only GLN184 is located in C-loop, 
and only LYS143 directly interacts with C-loop. Most of the 
residues are thus located apart from the region displaying the 
major drift during the transition. This is the sign of a domino effect 
relating the conformation of these residues to the C-loop opening 
[82]. 

The TAMD simulation produces spontaneously the holo-to-apo 
transition, accelerating by a factor of 1 0 the transition obtained as 
described above by altering the initial configuration in standard 
MD. The TAMD results give very qualitative estimation of tiie 
energy barriers which are crossed during the transition, through 
the observation of the effective energy values for which the 
conformational transition is observed. The value is in the range 3— 
5 kcal/mol for the AChBP transition, which is found to be in 
qualitative agreement with the free energy barrier of about 2 kcal/ 
mol between the open and closed C-loop states of the ECD of al 
nAChR, calculated by umbrella sampling using the LYS143/ 
TYR188 bond as collective variable [50]. In the presence of 
lobeline, the acceleration of the transition by TAMD is reduced, 
and the transition, not seen at ;8 '=3 kcal/mol, starts to be 
observed from equal to 5 kcal/mol. The increase of about 2 
kcal/mol induced by the presence of lobeline is in qualitative 
agreement with the affinity of 0.3 nM [27] experimentally 
measured for this ligand on Aplysia califomica AChBP. Interestingly, 
the TAMD acceleration also influences the early steps of lobeline 
dissociation. Indeed, the complete dissociation of lobeline from the 
binding pocket is achieved, for equal to 10 kcal/mol, in the 
nanosecond timescale. 

The analysis of the hydrogen bonds along standard MD and the 
TAMD trajectories permitted to describe the transition from the 
holo to the apo state in terms of the evolution of a network of 
interactions. Among the hydrogen bonds, the one between 
LYS143 and TRP188 has been shown previously in the literature 
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[50] to be important. But, also other hydrogen bonds like 
SER146/TYR93 have been found even more relevant to the 
mechanism of the transition, as pointed out by the TAMD 
simulations. Most of the residues involved in hydrogen bonds are 
conser\'ed in nAChRs and we can thus expect similar behavior for 
them in nAChRs. Their importance on the binding/gating events 
in nAChRs has been already tested by site direct mutagenesis 
[83,84], validating the model proposed here. 

All the observed conformational changes are quite local, but 
nevertheless relevant according to what it has proposed in the 
literatur(;, in particular by Auerbach and coworkers, in the so- 
called conformational wave model [85,86] where the ensemble of 
variations in local interactions is thought to produce the overall 
allosteric mechanism. The results obtained here on AChBP could 
be then fully relevant for the nAChRs. Furthermore, the effects of 
lobeline on tlie TAMD transition agrees with the lobeline function 
on nAChRs as fuU agonist, as the presence of lobehne decrease or 
suppress the acceleration of the C-loop opening. 

Materials and Methods 

Starting Structures 

It is known that Acetylcholine (ACh) is the natural ligand of 
AChBP, but, up to now, there is no X-ray crystallographic 
structure of AChBP bound to ACh. As it is quite important to start 
from a ligand position assessed by precise information from 
structural biology, a high-resolution structure of AChBP bound to 
lobeline was used, rather than a model obtained by docking ACh 
to an empty AChBP structure. Moreover, the choice of lobeUne as 
ligand is based on its therapeutic importance [8,9]. As starting 
conditions for the apo and holo AChBP conformation we used the 
X-ray crystallographic structures PDB entry names: 2W8E [33] 
and 2BYS [27] oi Aply.sia califomica AChBP. They were the highest 
resolution AChBP structures available (2.05 A for 2BYS and 1.9 A 
for 2W8E), at the time this work started (a new AChBP apo 
structure has been deposited later, PDB entry: 2Y7Y at 1.89 A). 
Missing residues in structures 2BYS (first loop after the helix al) 
and 2W8E (in loop C) were reconstructed using homology 
modeling with Modeller v9.8 [87,88]. In what follows, we refer 
to the original residue numbering in 2BYS [27]. Note that residue 
numbering is shifted by +2 with respect to the original numbering 
in 2W8E [33]. 

In order to investigate the role of selected residues on the stability 
of the AChBP conformations, we modified their arrangement in the 
initial protein conformations thus creating alternative starting points 
for the lobeline-bound AChBP structure. Details are fully given in 
Text SI (see Section "New AChBP models"). 

MD and TAMD Simulations 

MD simulations of AChBP were carried out for several systems in 
explicit water: the apo pentamer (P), the lobeline-bound pentamer 
in the presence (Pl+L) or in the absence (PI) of lobeline. The 
starting points for P and Pl+L/Pl simulations were respectively the 
first five chains, called A to E, of the PDB entries 2W8E and 2BYS, 
completed as explained. Starting from the new AChBP structures, 
additional MD simulations were carried out both in the presence 
and in the absence of lobeline. Set up of the systems and the 
standard MD simulation protocol are fully described in Text S 1 (see 
section "Molecular Dynamics simulations"). 

As for TAMD, the theoretical basis was originally presented by 
Maragliano and Vanden-Eijnden [58]. The coupled system of 
equations describing TAMD, for a system whose configuration is 
specified by xeR^, is: 



' ' (1) 

where /= 1,...3N and 7=l,...m. 6(x) = (6](x),...,6if,(x)) are 
collective variables (CVs) that are functions of the atomic 
Cartesian coordinates and z = {zi,...,z„) is a set of particular 
values of these variables; V{x) is the inter-atomic MD potential; k 
is the spring constant coupling the CVs with their target values; y 
and y are friction coefficient; 17 is a white-noise process, i.e. a 
Gaussian process with mean zero and covariance 

<^rii''{t)rij^(s) = SijS{t — s)y where a = x,z; j8= 77^ where ks is 

ksl 

the Boltzmann constant and the T is the system temperature; 

P= -. — where T is an artificial temperature, T>>T. The 
ksT 

aforementioned set of equations describe the motion of x{t) and 
9{t) over the following extended potential: 



7 m 

Uk{x,6)= {ej{x)-zj)\ (2) 

7=1 

The advantage of TAMD [58], is that if (i) y is chosen sufficientiy 
large so as to guarantee that z{t) evolves much more slowly relative 
to the fundamental variables, and (ii) k is sufiiciendy large such 
that Zy(;) « 6y(x(;)) at any given time, then the z variables 
effectively evolve according to the equation: 

Wjij = - ^ + ^jir'ymjnlit)- (3) 

where F{z) is the free energy associated to the CVs 9{x). 
Therefore the points in the CVs space are sampled according to 
the Boltzmann factor exp{ — PF) which is tuned to be more 
uniform than the associated physical Boltzmann factor exp{ — fSF) 
by taking T> >T. TAMD therefore accelerates the sampling of 
the conformational landscape through the CVs space by 
increasing the probability of visiting points with relatively low 
physical Boltzmann factors. 

The procedure for determining suitable values of y and spring 
constant k is described in Text SI (see Section "Determining the 
value of fictitious friction y and spring constant k", and Fig. S6). 
TAMD trajectories were initiated from: (i) the 40 ns frame of MD 
trajectories Pl-I-L and PI, (ii) the 40 ns frame of the Pl-I-L 
trajectory, to which the lobeline has been removed. 

As collective variables we used the Cartesian coordinates of the 
center of mass of the loop C (residues 183-194) and cys-loop 
(residues 128-140) of each subunit of AChBP, for a total of 30 
variables. We could do this because, when compared with other 
available enhanced sampUng methods, TAMD allows the use of a 
larger number of variables. Since we used a large number of 
collective variables and we observed the opening mechanism in 
the TAMD simulations, we did not compare alternative choices. A 
posteriori, the analysis of two different distances to assess C-loop 
opening, including the distance between the loop and the facing 
monomer surface (see below, "Trajectory analysis"), shows that 
although these distances are not directly evolved in TAMD, they 
sample the range of values from the closed to the open state. In 
particular, we have considered "open" those conformations that 
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showed a value of the distance similar to the one in the crystal 
structure 2W8E. 

TAMD simulations were run from 10 up to 30 ns, with 
temperature f of about 1500 K, 2500 K, 3500 K and 5000 K, 
corresponding to = 3, 5, 7 and 10 kcal/mol. In TAMD runs, 
the atomic force coming from the couphng potential and the 
evolution of the CVs are implemented via a Tel script linked to the 
NAMD code. 

A list of the MD simulations started from the original 
conformations and of the TAMD simulations is given in Table 1, 
together with other MD simulations started from modified 
conformations. Systems set-up details are given in Table S3. 

Trajectory Analysis 

In order to characterize the protein conformational changes, we 
used, following [89], three distances between the centers of mass of 
Ca atom groups. For the loop C, the distance d'"'™ (Fig. S5A) 
internal to each subunit, is calculated between the center of mass 
of residues 183-194 of loop C in the principal subunit and the 
center of mass of the residues 143 and 144. For the cys-loop, two 
distances d^' (Fig. S5B) and d°ff' (Fig. S5C) are defined, 
between the cys-loop center of mass (residues 128-140) and the 
center of mass of residues 179-181,198-204 and 120-127,49- 
60,30-43, found respectively in the inner and outer jS-sheet of the 
same subunit. These /J-sheet residues were picked up as the ones 
displaying average Mean Square Fluctuations smaller than 0.8 A 
in the P simulation (see Results section), to guarantee a "rigid" 
back-wall against which the cys-loop displaces. To study the role of 
specific residuc--r(;sidue interactions, an intra-protein hydrogen 
bond analysis has been performed along all MD trajectories (see 
Fig. 5, Fig. 7 and Fig. S7). Details are fuUy described in Text SI, 
see Section on "Hydrogen bond analysis". 

Supporting Information 

Figure SI Structure of lobeline. Schematic structure of the 

lobeline molecule. 

(TIFF) 

Figure S2 Root Mean Square Deviation of individual 

subunits. Root Mean Square Deviation (in A) of individual 
subunits calculated from the starting conformations (after 

/ ^ — ^ fef 2 

equilibration) as RMSD{tj) = \j ^^'-' — where XcStj) 

is the position of the atom at y'th time step, r'^i' is the position of 
the Ca atom in the reference structure, is the total number of 
Ca atoms in the subunit. Upper panel: P; central panel: PI; lower 
panel: Pl+L. The RMSD values are calculated after removing the 
roto-translational body motions of the single subunits [90] . The 
curves are colored according to the scheme in Fig. 1 in Main text. 
(TIFF) 

Figure S3 Distribution of d^"" values. Distribution of d^"" 
values (A) calculated over the TAMD trajectory at different time 
slices (see legends), at = 10 kcal/mol. Upper panel: (Pl)io; 
lower panel: (Pl-t-L)io. 
(TIFF) 

Figure S4 Protein residues involved in hydrogen bonds. 

Protein residues (in licorice) involved in hydrogen bonds analysed 
in this work. The protein is shown as a cartoon model. In panels C 
and D, the principal and complementary subunits are shown in 
gray and yellow, respectively. 
(TIFF) 



Figure S5 Schematic representation of the C-loop and 
cys-loop distances. Schematic representing of the C-loop and 
cys-loop distance parameters analyzed in this work. In A) the blue 
line indicates the d^""; B) d^f'; C) d^'f . Protein is shown as a 
cartoon model. 
(TIFF) 

Figure S6 Running average of the restraining force for 
each CV and 9j{x{ty) — Zj versus time. A) Running average of 
the restraining force for each CV, during MD simulation in which 

the CVs are fixed to their initial values; k= WOkcal/molA ^. B) 
6j(x{ty) — Zj values versus time for restrained dynamics,y= l,..,M 
where M is the total number of collective variables; 
k=mkcal/molA^. 
(TIFF) 

Figure S7 Hydrogen bond analysis. Number of subunits for 
which hydrogen bonds are present 50-90% of the time (A,C) or 
more than 90% of the time (B,D). The analyzed trajectories are 
MD (A,B) and TAMD at ^ = 3 or 5 kcal/mol (C,D). The residue 
pairs are: ASP159/ARG59 (black), GLU153/ARG79 (cyan), 
GLU193/ARG79 (magenta), GLU193/LYS25 (green). The 
hydrogen bond is considered to be formed when the minimum 
distance among all possible donor/ acceptor distances, for each 
pair, is less than 2.5 A. 
(TIFF) 

Figure S8 Probability of all atoms SER146 and TYR93 
RMSD. Probability of all atoms SER146 and TYR93 RMSD 
value less than 2 A from the reference states: P/1, P/2, Pl+L, P/1, 
& P/2, P/1 & P/2 & Pl+L, P/1 & Pl+L, P/2 & Pl+L, 
along the TAMD trajectories PI (black bar) and PI -40ns (dark 
gray bar) at ;8~ =3 kcal/mol. 
(TIFF) 

Figure S9 Probability of all atoms SER146 and TYR93 

RMSD. As in Fig. 19, along the TAMD trajectories (Pl)5 (black 

bar); (PI)? (dark gray bar) and (Pl)io (gray bar). 

(TIFF) 

Figure SIO Probability of all atoms SER146 and TYR93 
RMSD. As in Fig. 19, along the TAMD trajectories (Pl+L)3 
(black bar), (Pl+L)5 (dark gray bar), (Pl+L)? (gray bar) and (P1+ 
L)io (light gray bar). 
(TIFF) 

Table SI Distances between lobeline and AChBP atoms. 

(PDF) 

Table S2 Distances between lobeline and AChBP atoms. 

(PDF) 

Table S3 Size and contents of the simulated systems. 

(PDF) 

Text SI. 

(PDF) 
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